NYO-1480-69 


Courant  Institute  of 
Mathematical  Sciences 

AEC  Computing  and  Applied  Mathematics  Center 


A  Detached  Shock  Calculation  by 
Second  Order  Finite  Differences 

Arnold  Lapidus 


AEC  Research  and  Development  Report 

Mathematics 
February  1967 


New  York  University 


NEW  YORK  UNIVtRSITY 
COURANT  INSTJTUtE  -  LIBRARY 
2-5lMerJT^'     ^J.^Y..^(•  NY  i"f«ll 


UNCLASSIFIED 


AEC  Computing  and  Applied  Mathematics  Center 
Courant  Institute  of  Mathematical  Sciences 
New  York  University 


Mathematics  NYO-1480-69 

A  DETACHED  SHOCK  CALCULATION  BY  SECOND  ORDER 
FINITE  DIFFERENCES 

Arnold  Lapidus 


Contract  No.  AT(50-l)-1^8o 


UNCLASSIFIED 


New  YORK  UNIVERStTY 
COURANT  INSriTUTS  •  LflWART 


^ 


11 


Table  of  Contents 

Page 

I.  Introduction  1 

II.  Differential  Equations  4 

III.  Geometry  9 
rv.   Difference  Equations  15 

V.  Stability:   Simplified  Lax-Wendroff 

Artificial  Viscosity  52 

VI.  Computation  37 

A.  Test  Case  37 

B.  Initial  Flow  37 

C.  Verifications  39 

1.  Bernoulli  Steady  State  Constant  39 

2.  Stagnation  Point  Pressure  40 

3.  Standoff  Distance  of  the  Shock  as 

a  Function  of  Time  40 

k.      Pressure  Profile  on  the  Body  4l 

5.  Final  Shape  of  The  Shock  42 

6 .  Summary  42 

Graphs  44 

Bibliography  49 


Ill 


ABSTRACT 

A  detached  shock  problem  for  a  symmetric  curved 
convex  cylindrical  body  moving  parallel  to  its  plane  of 
symmetry  was  solved  using  a  third  order  accurate  Richtmyer 
form  of  the  Lax-Wendroff  conservation  equations.   One 
Innovation  is  an  easy  to  use  "artificial  viscosity"  term 
which  preserves  the  high  order  of  accuracy  of  the  calcula- 
tion while  removing  the  non-linear  instabilities  which 
otherwise  appear  in  the  shock  region  and  near  boundaries. 
Another  Innovation  is  a  simple  transformation  of  Cartesian 
space  which  changes  the  curved  body  into  a  straight  line, 
thus  reducing  the  large  number  of  special  points  and 
irregularly  shaped  mesh  regions  which  would  otherwise 
appear  in  the  difference  method  calculation.   Such  trans- 
formations are  shown  to  preserve  the  conservation  property 
of  the  system  of  differential  equations.   Other  aspects  of 
the  third  order  artificial  viscosity  term  and  the  trans- 
formation are  discussed.   The  results  of  a  numerical 
calculation  on  a  CDC  6600  Computer  are  compared  to  known 
results. 


I.   Introduction. 

Consider  a  smooth  plane-symmetric  convex  cylindrical 
body  moving  parallel  to  its  axis  of  symmetry  with  constant 
supersonic  speed  through  a  perfect  compressible  gas.   A 
steady  state  consists  of  a  detached  shock  at  some  distance 
from  the  body  which  has  a  shape  which  depends  on  the  shape 
of  the  body  and  its  speed.   There  is  a  region  of  flow  behind 
the  shock  and  in  front  of  the  body  in  which  the  flow  is 
subsonic,  and  a  region  behind  the  shock  in  which  the  flow 
is  supersonic,  these  regions  being  separated  by  the  so 
called  sonic  line.   The  problem  which  is  considered  here 
is  to  determine  the  position  of  the  shock  near  the  sonic 
region  and  the  flow  in  the  sonic  region  for  a  specific 
case  in  which  the  steady  state  configuration  looks  like 
Figure  1. 


shock 
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boundary  of  body 
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Figure  1 . 


The  detached  shock  problem  has  been  tackled  by 
diversified  and  Interesting  mathematical  methods. 
In  the  method  of  lines  of  Belotserkovskll  and 
Dorodnltsyn  [9],    the  region  of  Interest  Is  divided  Into 
strips  In  which  the  partial  differential  equations  are 
reduced  to  coupled  ordinary  differential  equations  that 
are  then  solved  numerically.   The  inverse  methods  ['i,?,!'^-] 
are  those  in  which  the  position  of  the  shock  is  assumed  to 
be  known  and  the  position  of  the  body  is  determined.   By 
iterating  the  shape  of  the  shock,  a  flow  can  be  found  in 
which  the  computed  shape  of  the  body  is  close  to  that  of 
the  problem.   Van  Dyke's  method  [l4]  is  very  fast  and 
accurate.   In  [?],  the  constructive  method  of  the 
Cauchy-Kowalewski  theorem  is  used  to  obtain  very  accurate 
numerical  results.   Garabedian's  interesting  method,  used 
in  [6],  is  to  continue  the  flow  into  the  complex  plane  in 
order  to  obtain  a  hyperbolic  initial  value  problem. 

The  detached  shock  problem  is  also  a  favorite 
practice  problem  for  developing,  testing  and  demonstrating 
time  dependent  methods  as  in  [5,8,15,16]. 

In  Bursteln's  work  [8]   a  time  dependent  flow  which 
tends  to  the  steady  state  is  calculated  using  conservation 
equations  and  the  shock  is  determined  as  a  region  of  rapid 
variation  of  the  flow  quantities;  we  proceed  similarly. 
The  new  features  here  compared  to  [8]  are: 


(1)  A  simpler  artificial  viscosity  term  Is  used. 

(2)  Since  the  boundary  of  the  body  Is  a  curve  Instead 

of  a  straight  line,  a  transformation  of  the  Cartesian 
plane  is  effected  which  maps  the  curved  body  onto  a 
straight  vertical  line,  and  the  finite  difference 
calculation  is  carried  out  in  this  non-physical 
coordinate  system.   The  conservation  laws  are 
rewritten  as  conservation  laws  in  this  new  system. 

(3)  To  continue  calculating  near  the  upper  boundary  of 
the  mesh,  flow  quantities  are  extrapolated  from  the 
interior  points  to  the  boundary. 

Finite  difference  methods  of  third  order  accuracy 
similar  to  those  in  [8]  were  employed  in  the  computation 
except  near  the  body  and  the  artificial  upper  boundary 
where  the  methods  used  were  only  second  order  to  keep  the 
computation  from  becoming  unstable.   The  stability  and 
accuracy  of  the  computation  will  be  discussed  in  the 
following  sections  on  the  differential  equations  as  well 
as  in  later  sections  on  the  numerical  results. 


II.  Differential  Equations . 

Fundamental  In  this  approach  to  solving  the  detached 
shock  problem  Is  the  use  of  the  conservation  form  [4,5] 
of  the  equations  of  time  dependent  two  dimensional 
compressible  fluid  dynamics: 
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are  cartesian  coordinates 

Is  time 

Is  mass  per  unit  volume 

Is  the  horizontal  velocity  component 

Is  the  vertical  velocity  component 

Is  p  •  u 

Is  p  •  V 

Is  total  energy  per  unit  volume 

Is  given  by  e  =  p/(r-l)  +  p{u^+v^)/2 

Is  the  ratio  of  specific  heats  (we  are  assuming 

a  gamma-law  gas ) 


In  this  approach  the  shock  is  not  a  boundary  but 
simply  an  Interior  region  of  rapid  variation  of  the  flow 
quantities.   In  fact  the  reason  for  using  the  conservation 
form  of  the  equations  is  that  it  Includes  the  shock 
conditions  as  well  as  the  differential  equations  of  motion 
of  the  gas.   Later  we  will  use  a  transformation  of  space. 

The  following  theorem  shows  that  if  a  non-singular 
transformation  of  space  is  made,  then  any  conservation 
law  expressed  in  the  old  coordinates  can  always  be  rewritten 
as  a  conservation  law  in  the  new  coordinates. 


Theorem  1.   Under  non-singular  space  transformations, 
conservation  laws  are  transformed  into  conservation  laws. 

Proof;   A  conservation  law  is  expressed  by 

(2)  Tu^dx  +   Tf^  dS  -  0 

G        dG 

where  P  =  F  .  n,  in  which  F  is  a  vector  (F, , . . -jF  ) 
and   n  =  (n,,...,ni  )  Is  the  normal  to  the  surface  G. 

After  a  change  in  variables  from  x, ,...,x  to 
i-]>  '  '  '  )ijj.   we  want  to  show  that  (2)  changes  into  an  equation 
of  the  same  form  in  the  new  coordinates.   Let 

J  =  S(x^,  .  .  .,Xjjj)  /  ^(^i»  •  • -^Ci^) 

and  let  J  be  the  Jacob  Ian  for  the  surface  element  of  (2) . 
s 


Let  M  be  defined  as  the  matrix  which  transforms  the 

direction  of  the  normal  v  to  the  surface  In  (f ,,..., ^  ) 

"1      m 

space  Into  the  direction  n,  of  the  associated  normal  to 
the  surface  In  (x.,,...,x^)  space.  The  Integral  (2) 
becomes 

(5)  /(Uj)d^  +  r(p.n)Jg  dS  =  0 

G         ^G 
where  bars  denote  the  quantities  In  the  new  coordinate 
system.   Now 


(p.n)Jg  =  (Jg  .  F,n)  =  (JgP,M-v)  =  (M'^JgF,v) 


where  M  means  transpose  of  M.   This  shews  that  (P  .  n)J 
is  the  component  of  a  vector  normal  to  'S  and  therefore 
(2)  and  (3)  are  of  the  same  form.   Q.E.D. 

A  practical  way  to  determine  the  components  of  the 
flux  is  given  by  the  following.   Use  the  differential 
form  (1)  of  the  conservation  laws.   Using  the  chain  rule 
and  multiplying  (1)  by  J  we  get 

m   m  Sf  .  hi, 

{5a)        (J  U)t  +  (  ^Z  ZZ  ^  ST  ^^  =  ^   ' 

^  j=l  k=l  ''^k  ''^J 

According  to  the  theorem,  the  spatial  part  of  (3a)  is  in 
conservation  form;  it  is  not  hard  to  show  that  it  is 


7 


(5b)         (J  U).  +  T"  (  V-  J  F.  -Jli^K   =0. 


(J  U).  +XI  (  IZ  J  F,  ^), 


It  Is  amusing  that  the  validity  of  (5b)  follows 
from  the  conceptual  argument  presented.   Of  course, 
(3b)  can  also  be  verified  by  a  computation.   We  note 
that  (5a)  can  be  written  as 


ju),  +^  ^^"i---^.i-i>^.i>^.m---\.^ 

But  according  to  a  theorem  of  advanced  calculus  [I3] 


S(x, ,  .  .  .,x._,  ,F.,x.  ,  ,  .  .  .,x  )  

6(i^,...     '^   '  '^  ;rT  ^^  °^  *^^  ^°^"^  ^  ^Pk^e  ^^^^  • 

Specifically,  the  transformation  which  we  used  was 

(4)  4^  =  x/h(y)    ^2  "  y  ^^  which  h'(y)  exists  and  h(y)?^0. 
This  changes  the  fundamental  equation  (1)  into 

(5)  (h(e2)U)^  +  (f-h'(e2)  •  eigi)|   +  (gh)^   --^0, 

and  maps  the  curve  x  =  h(y)  into  the  line  i^  =1 .      In  this 
problem  the  curve  x  =  h(y)  is  the  boundary  of  the  body. 
Such  a  transformation  of  the  space  of  independent 
variables  is  employed  as  a  means  of  simplifying  computation 
near  the  body.   It  is  obvious  from  Theorem  1  [eq.  (3b)]  that 
the  transformation  does  not  affect  the  existence  of  the 
shock  or  its  position. 
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In  summary,  we ' ve  shown  that  the  conservation  . 
properties  of  the  system  are  preserved  under  the  trans- 
formation and  that  the  essential  properties  of  the 
shock  are  also  preserved  under  the  transformation. 
It  Is  also  obvious  that  if  the  function  h  satisfies 
the  symmetry  property  h(y)  =  h(-y)  then  the  symmetry 
of  the  solution  to  the  transformed  system  is  also 
preserved. 


III.  Geometry. 

The  exterior  of  the  obstacle  in  which  the  flow  takes 
place  is  unbounded.   Since  only  a  finite  number  of  lattice 
points  can  be  employed  in  any  actual  computation,  the 
neighborhood  of  Infinity  has  to  be  dealt  with  in  some 
summary  fashion.   In  this  work,  we  resort  to  the  simple 
and  crude  expedient  of  confining  the  calculation  to  a 
finite  region  of  the  flow  field  which  contains  the  whole 
sonic  regime. 


Artificial  boundary 


boundary 

in   ^ J 

constant 
region 


consta 
region 


t 


line  of  symmetry 


boundary 
of  body 


Figure  2. 
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In  this  approach  the  region  Is  bounded  from  above 
by  an  artificial  boundary  line  on  which  boundary 
conditions  have  to  be  imposed  which  approximate  well  the 
state  of  affairs  along  this  line  in  the  true  solution. 
At  any  rate,  the  error  incurred  by  imposing  artificial 
boundary  conditions  should  be  comparable  to  discretiza- 
tion errors . 

The  remaining  geometrical  inconvenience  is  the 
Incommensurability  of  a  curved  body  with  a  rectangular 
mesh.   We  solved  that  problem  by  mapping  the  irregularly 
shaped  cut  out  configuration  onto  a  rectangular  region 
by  a  transformation  which  took  the  boundary  of  the  body 
into  the  right  side  of  a  rectangle.  The  mapping  which 
we  used  was 

|,  =  x/h(y)     ^2  =  y   where  x  =  h(y)  is  the  function 

which  traces  out  the  boundary  of  the  body  and  is  assumed 
to  be  continuously  differentiable,  symmetric  in  y  and 
non -vanishing.   Then  the  Jacobian  of  the  transformation 
Is  also  non-vanishing. 

Choose  the  coordinates  so  that  the  left  boundary  of 
Figure  2  is  the  line  x  =  0;  then  its  image  is  the. line 
^.,  =  0.  The  lines  y  =  0  and  y  =  .2  are  transformed  into 
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^P  =  0  and  ^2  =-2,  while  the  curve  x  =  h(y)  gets 
transformed  Into  the  line  ^,  =  1.   Figure  5  illustrates 
the  new  configuration  In  which  the  three  previously 
straight  boundaries  remain  so  while  the  curved  body  is 
transformed  into  a  line. 


artificial  boundary 

4-  >^ 


constant- 
boundary— > 


constant 


supersonic 


boundary 
of  body 


natural  boundary  of  symmetry 


Figure  3' 


This  transformation  also  saves  some  computer  space 
because  the  region  behind  the  body,  the  shaded  space  in 
Figure  2,  has  been  eliminated. 

The  mesh  which  was  chosen  for  the  machine  computation 
is  a  79  X  20  rectangular  grid  projected  onto  that 
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rectangular  portion  of  ^-j^^g  ^P^^e  given  by  [0,1]  ^  ^'^'Vj^p.yJ 
It  Is  the  set  of  points  (jA^^jkA^^^  ^^^^   J  =  0,...,78; 
k  =  0,  ...,19  and  M^  =   I/78,  Alj  =  ymax^^^*   ^^®  ^°'^^  ^° 
denoted  by  the  right  sided  boundary  of  the  mesh  since  the 
boundary  of  the  body  Is  transformed  Into  the  line  ^,  =  1. 
The  shock  Is  not  Introduced  as  an  explicit  boundary  but 
is  determined  as  the  locus  of  the  most  rapid  change  in 
the  values  of  appropriate  flow  quantities . 

In  summary,  a  region  of  interest  is  cut  out  of  the 
plane  and  then  transformed  into  a  rectangle.   The 
rectangular  region  is  then  represented  as  a  rectangular 
mesh  which  can  be  used  in  a  machine  computation. 
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IV.   Difference  Equations . 

In  this  section,  a  difference  scheme  for  the 
fundamental  equation  (l)  which  Is  similar  to  that  In 
[4]  Is  described  and  Is  shown  to  be  consistent  and  to 
have  truncation  error  O(a^)  In  the  smooth  part  of  the  flow 
which  Is  far  from  shocks,  body  or  other  boundaries.   We 
hoped  that  no  special  methods  would  have  to  be  used  at 
the  shock  except  to  use  the  two  step  scheme  and  that  the 
location  of  the  shock  could  be  determined  by  sharp 
variations  In  some  of  the  flow  quantities  such  as 
density,  entropy  or  pressure.   The  main  reason  for  this 
Is  that  the  two- step  scheme  Is  Itself  conservative  as 
shown  by  Lax  and  Wendroff  in  [5].   That  is,  solutions  U 
of  the  two  step  equations  satisfy  a  discrete  analogue 
of  the  integral  relation 

^^^         Iff  ^^t^  "^  w^P(U)  +  WyG(U))  dx  dy  dt  =  0  , 

where  w(x,y,t)  is  a  dlfferentiable  function  of  compact 
support.   Such  difference  schemes  are  conjectured  to 
satisfy  the  entropy  condition.   The  first  step  of  the 
two  step  scheme  is: 
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(7a) 


n+ 


" .  1,.,  i->5)c+":; 


J'   0>  '^'^      Q 


'"•3)c  +  ""rtk  +  "5k+i-^''5«k+i' 


At 


n 


n 


2Ax  ^^  1  ~  ^     1 ^ 

'^^   j+1  k+  I    J  k+  I 


At 


n 


n 


2Av  ^^    1     -  ^    1    ^ 

^^y    j+  A  k+i    j+  -|  k 


The  second  step  of  the  scheme  Is: 


Jk     Jk   ^   J+  ^  k   J-  "I  k 


n+1    n 
(7b)  U    =  U 


^y  J  k+  I  J  k-  -^ 


in  which 


(re)  u^^. 


K.1 

"jk 

"5. 

hJ 

f"  = 
Jk 


r  n 
'"jk 


/  n  \2  /  n  ,  n 
^'"jk)  /Pjk+Pjk 
„n  „n  /n 

L'^jk+Pjk'-^jk/Pjkj 


S. 


P"  - 

^jk- 


n   n  /  n 
"'jk'"jk/Pjk 
/„n  \2 /n    ,n 

/  n  ,  n  \„n  /  n 
>Jk+Pjk^"jk/Pjk 


and 


_n+  -i      n+  I 
F   ^  .  (F   2 


n+  3 
+  F   2    ^)/2. 

J+  "2  k-  "2 


(7d) 


_n+  "5      n+  -K        n+  "5 

G     1=^^   1    1+^   1    1^/2' 
J  k+  -^     J''"  "q  k+  "q    3~  "o         "o 


This  Is  a  nine  point  scheme  In  the  sense  that  only 
nine  points  of  the  mesh  at  time  n  'At  are  needed  to  obtain 
numerical  values  at  a  point  at  time  (n+1) At.   However,  four 
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additional  points  (A,B,C,D  of  Figure  h)    are  used  for  the 
computation  of  intermediate  values.   This  scheme  can  be 
used  only  for  points  which  have  all  eight  neighbors. 

Next  we  want  to  prove  consistency  for  this  scheme 
(equations  (7))  and  also  that  the  discretization  error 
Is  O(A^). 


(J-I,k-1) 


(J,k-1) 


(j+l,k-l) 


Figure  K 


Let  the  truncation  error  at  a  point  be  defined  as 
the  difference  between  a  computed  solution  and  the  true 
solution  having  the  same  initial  data.   This  is  also  called 
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discretization  error. 

Theorem  2.   The  scheme  described  by  formulas  (7)  Is 
consistent  with  the  fundamental  equation  (l)  and  has 
a  truncation  error  which  Is  0{^-^)    In  the  smooth  part 
of  the  flow. 


Proof:   Define  the  operators  N(U)  and  N.(U)  by  the 
equations 


N(U)  =  U^  +  Px  "^  ^y 

and  ,        , 

n+1   n  _^+"o    _i^+  o 

N.{U)  =  (U   -  U,.  )/At  +  (F   f   -  P   n       )/Ax  + 

^        Jk    ^^  J+^  k   J-  i  k 

_n+  "I     _n+  -5 

+  (G     1  -  G  i)/Ay  . 


In  order  to  prove  the  consistency  of  the  operators  N 

and  N.  It  will  be  shown  that 

n+l      n+l 
11m  |N   ^(U)  -  N.   "^(U)  I  =  0 
A-»-0  ^ 

for  arbitrary  smooth  functions  U.  First  we  assume  that 

11  1 

n+  -^     n+  -^  n+  -^ 

(8a)        (F   T   -  P   T   )/Ax  =  F^(U    )  +  0(A) 
J+-|  k   J-  f  k        ^  Jk 

and 

(8b)        (G   '^  T  G   ^  T)/Ay  =  G^(U..  '^)  +  0(A) 

J  k+  -i    J  k-  i  y  J^ 
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in  which  quantities  on  the  left  are  evaluated  by  the 
difference  scheme  and  the  quantities  on  the  right  are 
evaluated  by  the  assumed  known  function  U.   Then  we  have 

n+  ^  n+  ^        n+  I   n+  -J   n+  -i   u"^^-u" 

|N   '^  U  -  N.   ^  U|  =  |U,   ^+  F   '^+  G   "^ ■J^.,  ^^ 

'  A     '    '  t      X      y         At 

n+  3  n+  "5  n+  -i  n+  ^ 

-  F^  ^-G^      2-K)(A)|=  |U^   2-U^   2^(A)|  =  |0(A)| 

which  goes  to  0  as  A  -»-  0 . 

Now  we  still  have  to  prove  {8a)  and  (8b).  The  proof 
is  complicated  by  the  fact  that  the  approximations  to  the 
derivatives  in  the  difference  scheme  have  an  intermediate 
step.  Only  the  proof  of  (8a)  follows.  The  proof  of  (8b) 
would  follow  similar  lines.  To  make  the  proof  a  little 
clearer,  the  function  values  which  are  computed  from  the 
difference  scheme  are  denoted  by  Vs.   Then 

n+  —        n+  ^  n+  o 

(7d)         F   ^   =  (F(V   ^     )  +  F(v   (  .))/2 

J+-|  k       J+-^  k+^       j+-i  k--^ 

n+-5-  n    n 

(9)  F(V   I  .)    =   P(U   ,  U       U"    .  u" 

j  +  ^  k+^       Jk   J+-^'^   J  k+1    J+1  k+1 

where  "P  is  a  very  smooth  function  of  U.j^,U,^-j^  k' '  '  '  * 
We  are  assuming  that  U  is  a  smooth  function  of  x,y  and  t 

so  -j^ 

n+-p 

(10)  P(V   ^    -,  )  =  P(Ax,Ay)   where  P  is  a  smooth 

function  of  Ax  and  Ay 
and 
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.  2 
(11)  P(Ax,Ay)  =  P(0,0)  +  P^(0,0)Ax  +  P2(0,0)Ay+  P^^(0,0)-^ 

2 
+  P3_2(0,0)Ax  Ay  +  P22(0,0)-^  +  0{^^)     . 

Then  by  (7c3)  and  (11)  we  obtain 

_n+  ^   _n+  ^        n+  ^         n+  -J 
F   f   -F   f     F(V   f    t)+F(V   f    ,) 
J+-i  k   J-4  k      J4.  ^  k+ ^     j+  ^  k-  ^ 


Ax  5Sc 

(12)  ,  , 

n+  -i  n+  "5 

F(V   J    j^)+F(V   J    ^) 

-     ^ ^ ^ ^  =   2P^(0,0)  +  O(Ax^). 

Next,  we  will  show  that  2P,  (0,0)  =  F^(x,y, t )-fO(At) 
which  will  finish  the  proof  of  consistency.  From  (9) 
and  (7a) 

(13)  P(Ax,Ay)  =  F(|(U(x,y,t)  +U(x+Ax,y,t)  +U(x,y+Ay,t) 

+  U(x+Ax,y+Ay,t))  -  (At/2Ax) (F(U(x+Ax,y, t ) ) 
-  F(U(x,y,t)))-At/2Ay.  (G(U(x,y+Ay,  t))-G  (U(x,y,t))) 
Then  by  differentiating  (13)  with  respect  to  Ax  we  get 

(14)  P3^(Ax,0)  =   Fy- 

•[uy2-At/2Ax-    (FuU^-[F(U(x+Ax,y,t))-F(U(x,y,t))]/Ax)] 
and   for  Ax  -»-  0 
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P^(0,0)    =  Pg-[uy2   -   F^AtAl    =   P^(x,y,t)/2   +  0(At) 

This  completes  the  proof  of  consistency. 

Now  we  will  show  that  the  discretization  error  in 
the  smooth  part  of  the  flow  Is  0{^-^)  .      Let  the  solution 
of  the  differential  equation  be  U(x,y,t)  and  let  the 
result  of  calculating  the  approximation  to  the  solution 
at  time  t  +  At  from  values  of  the  solution  at  time  t  be 
denoted  by  VVi^  .   We  want  to  show 

(15)  U(x,y,t  +  At)  -  V^j^^  =  O(A^)  . 
From  the  differential  equation  we  get  that 

(16)  U(x,y,t+At)  =  U(x,y,t)  -  At(F^(x, t, t+At/2)+ 

+   G  (x,y,t+At/2))  +  O(A^) 


and  from  the  difference  scheme 


n+  i      n+  -J 


+1    TT/     J.^         At  rri"   2      h    2 

J+4  k    j-^  k 


(17)  V^J"-  =   U(x,y,t)  -H  Cp  I       -  F  (      )    - 


2 


J  k+-^    J  k-  -^ 

where  F  and  G  are  given  by  equation  (7d).   From  (l?)  It 
Is  clear  that  equation  (I5)  will  be  satisfied  if 

(18)   (P   ^   -  F   f   )/Ax  =  F^(x,y,t+At/2)  +  O(A^)  , 
J+-|  k    J— i  k 
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and  similarly  for  the  term  In  G(V).   To  prove  this, 
define  P(Ax,Ay)  so  that 


F(V"+^,     )  =  P(Ax,Ay)  = 
(U(x,y,t)-rtJ(x+Ax,y,t)-KJ(x,y+Ay,t)-KJ(x+Ax,y+Ay,t))A 


hdw          At/F(x+Ax,y,t)4-F(x+Ax,y+Ay,t)  F(x,y,  t)4-F(x,y+Ay,  t)  ^ 

^^^'  ^    "  2A5c^           2  2         ' 

At^G(x,y+Ay,t)-K}(x+Ax,y+Ay,t)  G(x,y,t)4G(x+Ax,y,  t)  ^ 

2Ay^           ^  2         / 


From  (19)  It  is  clear  that  the  equations 

n+  o  n+  p 

(20)  F{V   f    ,)  =  P(-Ax,Ay);  P(V   ^     ,)  =  P(A:c,-Ay) 
J-  -g  k+  -i  J+  -i  k-  -i 

n+-| 
P(V   J    ^)  =  P(-Ax,-Ay) 


hold.   If  we  assume  that  U{x,y,t)  is  sufficiently  smooth, 
we  can  expand  P  as 

(21)   P(Ax,Ay)  =  P(0,0)  +  P^{0,0)Ax  +  P2(0,0)Ay 

+  P^^(0,0)Ax^/2  4-  ?^^{0,0)^   Ay  + 
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From  (20)  and  (21)  It  then  follows  that 

n+  -5  n+  p 

(22)  [F(V    f     t)+P(V    f     n)/2Ax  - 
J+  i  k+  I      J+  -i  k-  ^ 

-  ^F(V    f     ,)  +F(V    f     t)1/2Ax]  =  2Pt(0,0)+O(A^) 
We  will  now  show  that 


(23)         Pi(0,0)  =  F^(x,y,t+At/2)/2  +  O(A^) 


and  then  (22)  implies  (18)  .   That  is,  it  is  sufficient 
to  show  (25)  in  order  to  show  that  the  discretization 
error  is  0{^^) . 

To  show  (25),  we  differentiate  (19)  with  respect  to 
Ax  and  obtain 


Pt(0,0)  =  lim  P„(U(x,y,t)-At/2  (F  -Kl  )  )  •  [U  (x,y,  t  )/2- 
•^        Ax-»^  0   ^  A  y     A 

Ay*  0 


Atrr.    t    ^A    ^'\^     At/F(x+Ax,y.t)-F(x,y,t)^   At  xy^ 
^^(F^(x+Ax,y,t;)+-5^( '—^  )-  -^  -^J 


Then  using  the  hypothesis  that  U(x,y,t)  is  a  solution 
to  the  fundamental  equation,  we  get 
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P^(0,0)=  Py(U{x,y,t+At/2))-  I  [U^(x,y,t)+At/2  (U^)^]+  O(A^) 
=  Py(U(x,y,t+At/2))/2.U^(x,y,t+At/2)  +  O(A^) 
=  F^(x,y,t+At/2)/2  +  O(A^).  '      Q.E.D. 


We  have  now  proven  that  this  nine  point  difference 
scheme  Is  consistent  with  the  fundamental  equation  (l) 
and  has  truncation  error  0(A-^^)  In  the  interior  parts 
of  the  mesh  which  are  not  in  the  shock  region.   The  left 
boundary  has  constant  values  and  those  at  the  axis  of 
symmetry  are  computed  by  reflection.   We  still  have  to 
describe  the  difference  schemes  which  are  applied  at 
the  two  non  trivial  boundaries:   the  body  and  the  upper 
artificial  boundary. 

At  the  upper  boundary  the  flow  quantities  were 
evaluated  by  linear  extrapolation  as  follows: 

n+1        n+1        n+1 

If  the  Mach  cone  at  every  point  of  the  upper 
boundary  were  to  leave  the  rectangular  region  of  calcula- 
tion, the  errors  introduced  by  extrapolation  {2k)    would 
not  affect  the  accuracy  of  the  computation.   While  the 
Mach  cones  do  not  leave  the  region,  the  curves  which  are 
the  envelopes  of  Mach  cones  starting  at  the  artificial 
boundary  end  at  points  on  the  body  which  are  still  in 
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the  supersonic  region.   The  numerical  experiments  show  that 
errors  which  are  Introduced  at  the  upper  boundary  affect 
the  subsonic  region  only  slightly.   The  linear  extrapolation 
(24)  Is  In  a  45  degree  direction  which  was  roughly  the 
same  as  that  of  the  nearby  characteristics.   Scheme  (24) 
was  unstable  at  the  early  stages  of  calculation  and 
therefore  was  replaced  by  (25)  where  the  factor  r  was  only 
gradually  built  up  to  1. 


(25) 


U    =  U        "^  ^*^"l-l  k-1"  " 
^Jk    ^J-1  k-1      J  -L  k  ± 


) 


J-2  k-2 


In  which 
r  =  .5   if  the  number  of  cycles  was  less  than  100  and  >  0 


r  =  .6 
r  =  .7 
r  =  .8 

r  =  .9 
r  =1.C 


It   II 


200  "  "99 

500  "  "199 

400  "  "299 

500  "  "399 


OD   "  "499 


The  top  point  of  the  body  is  also  treated  In  this 
ad  hoc  fashion  of  a  point  on  the  upper  boundary  but  the 
other  points  of  the  body  are  treated  in  a  manner  which 
Is  very  much  In  the  spirit  of  the  conservation  law 
approach.   To  explain  the  scheme  used  at  the  body  we 
will  refer  to  Figure  5. 


2k 


3  0 


k  0 


5  e 


^ 


,10 


Of-- 


d)' 


^ 


^o? 


1/  2 


ii  1 


1-^ 


<D"6 


boundary  of  body 


G  denotes  mesh  point 

A  denotes  Intermediate  point 


Figure  5.   Difference  Scheme  at  the  Body. 


The  flow  quantities  are  assumed  to  be  known  at  the 
points  labelled  1,2,3,4,5  and  6  at  time  t.   The  problem 
Is  to  obtain  values  of  the  flow  quantities  for  the  point 
labelled  1  at  time  t  +  At . 

We  think  of  the  flow  quantities  at  point  1  as  being 
affected  only  by  the  fluxes  P  and  G  across  the  dotted 
lines.   To  obtain  good  estimates  of  F  and  G  across  those 
lines  we  first  make  estimates  of  the  flow  quantities  at 
the  points  labelled  7,8,9,10  at  time  t+At/2. 
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At  7  (and  10)  we  calculate  as  at  ordinary 
intermediate  points.   That  is,  for  the  point  7  we 
consider  the  rectangle  whose  vertices  are  the  points 
1>^^5^6.   To  compute  the  change  in  the  flow  quantities 
at  point  7,  we  compute  the  sum  of  the  fluxes  through  the 
sides  of  the  rectangle  and  add  it  to  the  flow  quantities. 
At  time  t,  the  average  of  the  flow  quantities  at  points 
1,4,5  and  6  was  used  as  an  estimate  for  the  values  at 
point  7- 

A  point  such  as  8  (or  9)  has  to  be  treated  in  a 
special  way.   The  method  chosen  is  to  estimate  its  flow 
quantities  at  time  t  by  averaging  those  at  points  1  and  6. 
Since  point  8  is  in  the  same  large  rectangle  (1,4,5,6)  as 
point  7,  the  change  in  the  flow  quantities  due  to  fluxes 
Is  taken  to  be  the  same  as  at  7- 

Having  computed  the  flow  quantities  at  time  t+At/2 
for  the  points  7*8,9>10  we  use  these  values  to  compute 
the  flux  across  the  dotted  lines  to  give  the  change  in 
the  value  of  the  flow  quantities  at  point  1  between  time 
t  and  t+At.   The  formulas  which  follow  express  the  various 
quantities  involved  in  the  procedure  described  above. 


(26a)  Flux  into  rectangle  (1,4,5,6)  =  B„q 

=  -  ^   ((^6+  P?)/2  -  (F!;+F^)/2)-T^(y) 
'   -  2^  ((^6-^1^/2  -  (Gi;-KJ^)/2)-T2(x,y) 
-  2§  (^G;4g5')/2  -  (G^-K}^)/2) 
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n+  ^   -,    „   „   „   „     n+  -| 


(26b)        v^  2  =  I  (U^l+Uj+U^+Ug)  +  B^g 
(26c)        Vq  2  ^  1  (yn  ^  yn^  ^  B^^  2 


where  T^(y)  =  l/h(y)   T2(x,y)  =  -xh ' (y)/h(y ) . 

The  same  formulas  with  the  appropriate  changes  in 
subscripts  (referring  to  points  in  the  figure)  are  used 
to  obtain  intermediate  values  at  the  points  numbered  9 
and  10.   At  the  point  numbered  1  the  difference  formula 
Is: 

n+1    n   2At    ""*""!    ^"*'  "I        n+  ^  n+  ^ 


U,   =  U,  - 


Ax 


((F^   ^  +  Fg   ^)/2  -  (F^   ^+F^Q  '')/2).T^(y) 


11  11 

(27)       -T2(x,y)(^)((G^  2+G8  2)/2-(G^  ^+G^q^)/2) 


-{^){{G^+G^q)/2   -    (Gg+G^)/2)  . 


This  difference  scheme  is  consistent  but  has  an  error 
of  second  order  rather  than  a  third  order  error  as  at  ordinary 
points . 

The  physical  boundary  condition  is  that  the  flow  should 
be  parallel  to  the  wall  at  body  points.   This  was  not  used 
so  far;  we  impose  it  now  by  replacing  the  momentum  vector 
calculated  at  boundary  points  by  one  which  has  the  same 
magnitude  and  the  right  direction.   Denoting  by  m'  and  n' 
the  new  momenta,  we  have 


27 


(28)  a.  =    V(m"+^)'^  4-  (n"^^''^ 

where  (^,ti)  are  the  components  of  the  normal  to  the  body 
at  the  point. 

At  the  nose  of  the  body,  the  syinmetry  of  the 
configuration  is  used  to  reduce  the  point  to  an  ordinary 
point  of  the  body.   This  is  the  same  procedure  used  at 
other  points  of  the  axis  of  symmetry,  i.e.  the  point  is 
reduced  to  a  known  type  by  continuing  the  flow  so  that 
the  point  gets  its  full  complement  of  neighbors. 

We  have  now  explained  the  difference  method  used  in 
the  mesh  and  at  boundaries  and  also  discussed  the  consistency 
truncation  error  of  the  equations.   The  transformed  equations 
were  used  in  the  conservation  form  as  well  as  in  the 
(analytically)  equivalent  form: 


(29)         U^  +  F^-h'^(y)  -  x-G^-h'(y)-h"^(y)  +  Gy  =  0 


A  natural  question  to  ask  is  whether  the  divergence 
free  conservation  form  of  the  equations  would  be  more 
suitable  for  computing  than  the  equivalent  form  (29).   To 
Investigate  this  question  the  now  classical  equation 

(30)  U^  +  (u2/2)^  =  0 

was  used.   The  problem  chosen  for  experimentation  was 
to  solve 
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(31)    U^  +  (U^/2)^  =  0   on  X  e  [-Tr/k,TT/k) ,    t   e  [0,tt/2], 
with  initial  conditions 


(32)  U(x,0)  =1  for  X  <  0 

U(x,0)  =0  for  X  >  0. 

The  solution  to    (31)   and    (32)    is 

(33)  U(x,t)    =1      if     X   <    t/2,      0    <    t    <  7r/2; 

U(x,t)    =0      if     X  >    t/2    . 

We  now  consider  a  change  in  the  space  variable 
given  by 

(34)  i   =   sin  X,  -ir/h   <  x  <  v/k,    {-  Y^/2<i<'/^/2)  . 

The  transformed  function  U(^,t)  is  then  to  have  initial 
conditions 

U(^,0)  =  1  for  ^  <  0 

(35) 

U(^,0)  =0  for  ^  >  0 

and  it  will  then  satisfy 

^t  +  f^'/2)x  ^x  =   ° 
where 

i     =  cos  X  =  Vl  -  sin  X  =  yl-i 
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Dropping  the  bar  and  denoting  (1-^  y-f  ^  ^y   rp^ 
the  problem  In  the  new  variable  becomes 

(36)  U^  +  (U^/2)^  T  =  0 
with 

U(e,0)  =1   for  i    <  0 
and 

U(4,0)  =  0   for  ^  >  0. 

The  solution  will  then  be 

U(4,t)  =  1  for  arc  3ln  i    <   t/2, 
U(4,t)  =  0   for  arc  sin  4  >   t/2. 

Now  let  U/T  =  V;  then  the  problem  becomes 

(37)  V^  +  (tV/2)^  =  0 

with 

V(^,0)  =  1/T  for  i    <  0 
and 

Y{i,0)   =0    for  I  >  0. 

The  solution  Is 

V(e,t)  =  1/t  for  arc  sin  ^  <  t/2, 

V(^,t)  =  0    for  arc  sin  i   >   t/2. 
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Equation  (57)  Is  In  conservation  form  while  (56)  is 
not.   The  pertinent  question  is  whether  there  Is  any 
difference  in  the  results  of  using  a  reasonable  difference 
approximation  for  these. 

The  difference  equations  used  where 

n+1  n  n  -,.,  n2n2 

(38)  Uj        =    (Uj^j   +  Uj.j)/2   -    If    {(tlj^^)    -    (Uj.j)    )T(5) 

n+1  n  n  ,.,  n22  n22 

(39)  Vj        =    (Vj^,    +  Vj.,)/2   -  I  f    ((Vj^,)    Tj^,-(Vj_,)    Tj_,) 


The  results  showed  that  (38)  gave  five  place  accuracy 
except  in  the  shock  region  while  (59)  gave  only  three  place 
accuracy.   The  reason  for  this  difference  is  that  U  is 
constant  while  V  is  not,  therefore  the  error  committed  by 
taking  the  average  in  the  first  term  on  the  right  is 
negligible  in  (58)  and  appreciable  in  (59)- 

Schemes  (58)  and  (59)  both  had  second  order  accuracy. 
For  third  order  accurate  schemes  analogous  to  the  two  step 
method  in  two  dimensions,  computational  experiments  showed 
almost  no  difference  in  the  numerical  results-   The  differ- 
ence analogue  of  (29)  was  used  as  well  as  the  analogue  of  (5) 
for  the  fluid  dynamics  equations  in  two  dimensions.   The 
results  Indicate  that  (5)  gives  much  better  accuracy  for  these 
equations  than  (29)  does. 

In  uhls  chapter  we  discussed  the  difference  equations 
used  far  from  boundaries  and  their  consistency  and 
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truncation  errors.   We  also  discussed  the  methods  at 
the  upper  artificial  boundary  and  at  the  body.   The 
results  of  a  computational  experiment  for  a  one  dimensional 
problem  Justified  trying  a  slightly  different  form  of  the 
difference  equations.   Stability  will  be  discussed  In  the 
next  chapter  and  results  of  a  machine  computation  will  be 
analyzed  in  the  last  chapter. 
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V.   Stability:   Simplified  Lax-Wendroff  Artificial  Viscosity, 

The  major  difficulty  in  this  problem  was  maintaining 
numerical  stability.   At  various  times  instabilities 
appeared  near  the  body,  near  the  shock  and  at  the  upper 
boundary.   To  counteract  this  instability  we  introduced  a 
new  type  of  stabilizing  transformation  which  consisted  of 
replacing  the  flow  quantities  U    calculated  in  the  last 
section  by  new  ones  U'   obtained  by  smoothing  first  in 
the  i2_'^±rectlon   and  then  in  the  ^pdlrection  according  to 
the  following  prescription. 

(40a)     Ujfl  =  Uj^'  +  AC  A'dA'u^;^  kl'^'^Vi  ^)) 

(40b)     Uj-^  =  U--^  .  .0  .  A"(  I  A"  V  .:^\|  .A"  (U -^\)  ) 

Where  A'U^j^  =  U^^-U^.^  ^;  A"  U^^  =  U^^-  Uj^_^ 
and  C  is  a  constant  (taken  as  4 .  in  actual  computations) 
while  u  and  v'  are  horizontal  and  vertical  fluid  velocity 
components . 

Equation  (40a,b)  are  fractional  steps  [20]  for  the 
numerical  solution  of  the  diffusion  equation 


U^  =  ac(a5)[(|u^|u^)^  +  {\V^  |U^)^] 


Prom  this  equation  it  is  clear  that  smoothing  is  of  third 
order  and  consequently  does  not  affect  the  truncation  error 
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of  the  difference  scheme. 

Adding  such  a  term  to  the  equations  of  motion  is 
similar  to  the  von  Neumann-Rlchtmyer  approach  used  in  [l8] 
and  was  motivated  by  the  higher  order  artificial  viscosity 
of  Lax  and  Wendroff  [5]. 

In  [5]  Lax  and  Wendroff  introduced  a  class  of 
difference  approximations  to  differential  conservation 
laws  which  themselves  are  in  discrete  conservation  form. 
We  now  describe  this  class  and  for  the  sake  of  simplicity 
we  do  it  in  only  one  space  variable. 

Let  the  differential  conservation  law  be 

^t  "  -^x     where  P  =  P(U)  , 

with  Initial  conditions  U(x,0)  =  (|)(x).  Let  g  be  a 
function  of  21   arguments  which  has  the  property  that 

(41)  g(U, ...,U)  =  P(U). 

Consider  the  following  difference  approximation  to  the 
equation  U.  =  P  : 

(42)  Av/At  =  Ag/Ax 

where 

(43)  Av  =  v(x,  t+At)  -  v(x,t)  ,    v(x,0)  =  ^x) 
and 
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(44)  Ag  =  g(x  +  Ax/2)  -  g(x-Ax/2) 
with 

(45)  g{x+Ax/2)  =  g(v^+i^  v^+2*  •  •  ■''^^^ 

where  the   v's    are  at  points  distributed  symmetrically 
around  x  +  Zix/2 . 

Multiply  equation  (42)  by  a  smooth  test  vector  of 
compact  support,  integrate  with  respect  to  x  and  sum  over 
all  values  of  t  which  are  integer  multiples  of  At.   On  the 
left  side  apply  summation  by  parts,  on  the  right  side 
replace  the  variables  of  integration  by  x+  Ax/2  and 
x-Ax/2  in  the  two  Integrals  which  appear  there.   If  g(x) 
is  defined  at  non-mesh  points  by  linear  interpolation  from 
g(x+Ax/2)  and  g(x-Ax/2);  we  get: 

(46)  _^y^v;(x,t)-w{x,t-At)  ^(^^t)dx  At  -  jT  w(x,0)(l,(x)  dx 

=  -  JZ  f  ^d2^±^^l^^Li2^iMM.   g(x)  dx  At  . 

This  leads  to  the  following  Lax-Wendroff 

Consistency  Theorem.   Let  v  be  the  solution  of  the 
difference  equation  (42)  in  conservation  form  with  initial 
conditions  ^   and  suppose  that  as  At, Ax  ->  0   v  tends 
boundedly  almost  everywhere  to  some  limit  U.  Then  U  is  a 
solution  of  the  differential  conservation  law  and  has 
Initial  values  ({>. 
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This  consistency  theorem  Is  actually  central  to  the 
way  that  we  are  computing  the  shocked  flow.  We  already 
know  that  solutions  to  the  Integral  equation 

(47)        ff  ^^t^-^x^''   ^^  ^^   +  rw(x,0)(f){x)  dx  =  0 

satisfy  the  differential  equations,  the  Initial  conditions, 
and  the  shock  conditions.   The  consistency  theorem  says 
that  solutions  to  (42)  approximate  the  solutions  to  (47) 
near  the  shock  as  well  as  In  the  smooth  region.   Then, 
to  compute  shocked  flows  we  merely  have  to  Iterate  flow 
quantities  using  an  equation  of  the  class  defined  by  (42). 

We  show  now  that  If  we  apply  smoothing  to  solutions 
of  difference  equations  In  conservation  form,  the  smoothed 
functions  also  satisfy  difference  equations  in  conservation 
form.  More  precisely: 

Smoothing  Theorem.  If  a  difference  scheme  satisfies  a 
conservation  law,  and  if  smoothing  is  done  by  adding  a 
term  which  is  a  first  difference  of  a  function  S  which 
has  the  two  properties: 

1)  S  =  s(u_^,u_^^3^,...,u^_^) 

2)  S(U,  U,  . ..,  U)  =  0 

Then  the  scheme  with  smoothing  also  satisfies  the  same 
conservation  law. 
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The  proof  follows  simply  by  constructing 

and   using  the  difference  scheme: 

Av  ^  ^1 
At    Ax 

It  Is  clear  from  the  definition  of  S  that 
g- (U,U, . . . ,U)  =  F(U)  Implying  that  the  consistency 
theorem  holds  for  g,  as  well  as  for  g. 

Such  smoothing  has  been  used  by  Kasahara  [21]  In 
atmospheric  fluid  dynamics  problems  and  by  Rusanov  [22] 
In  a  variety  of  time-dependent  problems.  Kasahara 
reported  that  he  used  second  order  smoothing  once 
every  forty  cycles  to  maintain  stability  Instead  of 
third  order  smoothing  at  each  step  as  In  this  report 
We  tried  a  calculation  in  which  we  smoothed  at  alternate 
cycles.   It  becajne  unstable  very  quickly. 

In  this  section  we  have  explained  the  reasons  for 
using  smoothing.   We  have  shown  that  the  smoothing  retains 
all  the  accuracy  that  the  original  system  had  and  that 
consistency  in  the  sense  of  the  Integral  operator  is  also 
preserved.   In  the  next  section  the  results  of  a  machine 
computation  at  Mach  6.  will  be  discussed. 
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VI.   Computation. 

A.   Test  Case. 

The  flow  chosen  for  a  test  case  was  one  computed 
by  Eva  Swenson  and  reported  by  her  in  [6].   This  flow 
had  plane  symmetry  and  included  a  detached  shock  whose 
shape  was  prescribed  to  be 


X  =  3  /I  +  y 


/T^ 


where  x  is  the  distance  along  the  axis  of  symmetry  and 
y  is  the  distance  from  the  axis  of  symmetry. 

The  gas  was  assumed  to  be  perfect  with  gamma  =  1.4. 
The  Mach  number  at  od  was  6.   This  inverse  problem  was 
solved  to  five  figures  of  accuracy  by  Garabadian's  method 
of  characteristics  in  complex  space  [19]»  furnishing  the 
body  used  in  our  calculation. 

B.  Initial  Flow. 

Our  time  dependent  method  requires  initial  values 
of  the  flow  quantities  at  interior  points  as  well  as  at 
the  body  points.   The  choice  of  this  initial  flow  is  to 
a  certain  extent  arbitrary.   Presumably,  the  flow  will 
eventually  settle  down  to  the  correct  steady  state  no 
matter  what  data  one  starts  with.   Of  course,  the  closer 


Since  the  position  of  the  body  was  not  known  at  precisely 
the  points  needed  by  our  finite  difference  method,  we 
used  quadratic  interpolation  to  obtain  the  needed  values. 
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our  Initial  guess  comes  to  representing  the  stationary 
flow,    the  faster  It  will  converge    to  a  steady  state. 
In  our  calculation  the  Initial  data  at  the  body  were 
the  results  computed  by  Swenson;  those  at  the  upstream 
boundary  were  assigned  the  data  at  oo  .   The  Initial  flow 
quantities  at  all  other  points  were  arbitrarily  set  to  be 
linear  functions  of  i  ,    having  the  above  data  at  the  body 
and  at  the  left  boundary. 

It  was  noticed  that  with  initial  data  in  which  the 
pressure  and  density  at  the  body  were  set  equal  to  their 
values  at  oo  ,  the  calculational  scheme  became  unstable 
after  a  few  cycles.   A  possible  explanation  for  this 
instability  is  that  these  initial  values  are  too  far 
from  the  steady  state  solution  and  initially  give  rise 
to  too  violent  a  flow.  A  similar  difficulty  was 
encountered  by  Bursteln  [11].  A  way  out  of  the  difficulty 
was  to  iterate  a  number  of  cycles  using  only  the  first  step 
of  the  two  step  method.   This  procedure  is  known  to  be 
stable.   These  nunerlcal  results  Indicate  that  if  fluxes 
are  too  large,  the  two  step  method  is  unstable.   This  is 
an  instance  of  nonlinear  instability  because  for  linear 
systems,  stability  is  not  influenced  by  the  size  of  the 
initial  data.   At  present,  I  know  of  no  explanation  for 
this  instability  nor  for  the  subsequent  stability  of  the 
scheme  when  applied  to  the  more  accurate  data. 
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The  usual  stability  analysis  based  on  v,  Neumann 
criteria  is  not  applicable  here  because  in  that  theory 
one  has  to  assume  that  flow  quantities  at  neighboring 
lattice  points  differ  by  0(h),  whereas  in  our  calculation 
this  instability  occurred  in  the  shock  region  where  the 
flow  quantities  vary  rapidly. 

C .  Verifications . 

In  this  section  we  compare  the  results  of  our 
calculation  with  the  exact  steady  state  known  from 
E.  Swenson's  work.  We  also  present  various  flow 
quantities  as  functions  of  time. 

1.  Bernoulli  Steady  State  Constant. 

In  [1,  p.  300]  it  is  proven  that  for  a  steady  flow 
the  quantity 

B  =  (u^+v^)/2  +  1/(7-1) -C^ 

is  a  constant  along  every  streamline  even  when  the 
streamline  crosses  a  shock.   Since  each  streamline  comes 
from  upstream,  B  is  a  constant  throughout.   The  initial 
average  of  B  is  29-04.   The  final  average  of  B  is  28.55. 
The  correct  upstream  value  is  28.7OO.  The  average  of  B 
as  a  function  of  cycle  number  is  plotted  on  graph  Ql. 
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2.  Stagnation  Point  Pressure. 

The  graph  G2  Illustrates  the  pressure  at  the 
stagnation  point  as  a  function  of  cycle  number.   The 
pressure  at  cycle  number  0  Is  the  steady  state  theoretical 
pressure  p  consistent  with  the  parameters  of  the  problem. 
Referring  to  the  graph.  It  Is  seen  that  the  pressure, 
which  starts  out  at  p  ,  drops  very  rapidly,  then  turns 
around,  overshoots  p  and  finally  approaches  It 
asymptotically. 

3.  Standoff  Distance  of  the  Shock  as  a  Function  of  Time. 

The  equations  of  motion  assume  that  the  gas  appearing 
In  the  computation  is  a  "gamma  law"  gas,  that  is  that 

P  =  Ap^ 

where  A  Is  a  function  of  entropy  S  alone.   Specifically 
[1,  p.  10] 

A  =  ke^   (  =  1.   upstream) 

where  k  and  g  are  positive  constants.   Thus  A  is  a  monotonia 
Increasing  function  of  S  and  can  be  used  as  a  measure  of 
the  entropy  Jump  across  the  shock.   Using  this  quantity, 
a  simple  method  for  determining  the  position  of  the  shock 
would  be  to  march  from  in  front  of  the  shock  (the  constant 
region)  along  a  horizontal  line  until  a  point  is  reached 
at  which  A  is  greater  than  1  and  to  assign  the  position  of 
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the  shock  to  that  point.   Unfortunately,  due  to  oscillation 
of  the  flow  quantities,  the  quantity  A  oscillates  below  1, 
then  above;  then  sometimes  the  cycle  repeats.   Because  of 
these  Irregularities,  the  position  of  the  shock  is  taken 
to  be  halfway  between  the  last  point  of  an  unbroken 
sequence  of  points  at  which  A  =  1  and  the  first  point  of 
an  unbroken  sequence  of  points  with  A  >  1.   The  position 
of  the  shock  on  the  line  of  symmetry  is  graphed  in  G5  using 
this  method,  along  with  its  theoretical  position  in  units 
of  meshpolnts  from  the  body.  The  position  of  the  shock 
starts  at  61.5,  then  it  moves  lower  and  finally  returns 
to  a  position  around  57.5-  The  true  position  is  44  In  these 
units.   The  body  is  at  0. 

4.   Pressiire  Profile  on  the  Body. 

In  graph  G4,  the  pressure  profile  along  the  body  is 
Illustrated  at  three  different  times  and  the  pressure 
curve  obtained  by  Interpolation  from  [6]  is  plotted  on 
the  same  graph.   As  can  be  seen  there,  the  pressure  profile 
on  the  body  at  100  cycles  moves  far  below  its  exact 
position;  then  at  5OO  cycles  far  above,  then  comes  much 
closer  to  theoretical  results  at  1000  cycles.   The 
pressure  near  the  nose  of  the  body  (lower  point  number) 
is  relatively  more  accurate  than  it  is  downstream. 
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5.  Final  Shape  of  the  Shock. 

In  graph  G5,  the  position  of  the  shock  at  1000  cycles 
is  illustrated  and  compared  with  its  exact  position  and  the 
position  of  the  shock  obtained  using  a  finite  difference 
analogue  of  equation  (29).   It  is  clear  that  the  conservation 
analogue  gives  much  better  results  than  equation  (29) 
especially  near  the  upper  boundary. 

6.  Summary. 

In  summary,  the  shape  of  the  shock  is  close  to  the  exact 
shape  and  the  computed  stagnation  pressure  is  very  close  to 
its  theoretical  value.   The  shape  of  the  computed  pressure 
profile  is  close  to  the  shape  of  the  theoretical  one  although 
the  values  of  pressure  are  up  to  30°/)  in  error.   These  results 
are  much  better  than  those  obtained  from  the  non-conservative 
calculation  in  which  the  position  of  the  shock  deteriorated 
as  h'(y)  became  large  and  in  which  the  shape  of  the  pressure 
profile  differed  drastically  from  its  theoretical  shape. 

In  the  future,  we  will  make  use  of  conservation  laws  and 
their  conservative  analogues  as  well  as  artificial  viscosity 
and  space  transformations  in  truly  non-stationary  shock  wave 
problems. 
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